using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.

MANOVA for species level

I cannot do site contrasts because of our uneven samplign scheme. The MANOVA throws out samples.

fit_full_species_man$pca_summary 
## Importance of first k=7 (out of 35) components:
##                           PC1    PC2     PC3     PC4     PC5     PC6    PC7
## Standard deviation     1.0736 0.7112 0.47335 0.37674 0.33980 0.30966 0.2955
## Proportion of Variance 0.4066 0.1784 0.07903 0.05006 0.04073 0.03382 0.0308
## Cumulative Proportion  0.4066 0.5850 0.66402 0.71409 0.75481 0.78864 0.8194
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex, data = data), test = "Wilks")
##               Df   Wilks approx F num Df den Df    Pr(>F)    
## hostRace       4 0.03980   40.088     28 776.62 < 2.2e-16 ***
## sex            1 0.81535    6.956      7 215.00 1.795e-07 ***
## hostRace:sex   4 0.79391    1.833     28 776.62  0.005633 ** 
## Residuals    221                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PCA for Species Level

LDA for Species level

##            Reference
## Prediction  Cingulata Cornivora Mendax pom Zepheria
##   Cingulata        25         0      0   0        0
##   Cornivora         1         5      0   0        0
##   Mendax            0         0     60   5        1
##   pom               0         0      1 116        0
##   Zepheria          0         0      0   1       16
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   9.610390e-01   9.390180e-01   9.273307e-01   9.820323e-01   5.281385e-01 
## AccuracyPValue  McnemarPValue 
##   1.508463e-49            NaN

Bayesian Model

Collapsed all 35 chc measurements into a value column with a corresponding chc column detailing what chc was measured. Fitted a model to evaluate the efect of host and sex on chc values. The model accounts for uneven site effects and individual effects.

All Bayesian models were created in Stan computational framework (http://mc-stan.org/) accessed with brms package (Bürkner, 2017). To improve convergence and guard against overfitting, we specified mildly informative conservative priors. We ran 4 chains with 2000 interations. We checked for convergence and found all models mixed well. (Pareto K estimates are all good!)

#1 value ~ hostRace * sex + (1|site) +(1|chc) + (1|ID)

model comparison?

Model Converged

mod = readRDS("fit_chc_model1.rds")
plot(mod)

model fit is good = posterior values match real values

brms::pp_check(mod, resp ="value") #great
## Using 10 posterior samples for ppc type 'dens_overlay' by default.

Model fit is good = r2 == 88%

bayes_R2(mod, summary = TRUE)
##     Estimate    Est.Error      Q2.5     Q97.5
## R2 0.8790622 0.0009515364 0.8771306 0.8808946

model summary

tab_model(mod)
  value
Predictors Estimates CI (95%)
Intercept -0.25 -0.49 – -0.00
hostRace: Cornivora -0.13 -0.29 – 0.03
hostRace: Mendax 0.17 0.09 – 0.26
hostRace: Zepheria 0.33 0.23 – 0.43
hostRace: pom 0.11 0.04 – 0.17
sex: Male 0.03 -0.06 – 0.11
hostRaceCornivora.sexMale -0.14 -0.35 – 0.06
hostRaceMendax.sexMale -0.07 -0.17 – 0.03
hostRaceZepheria.sexMale -0.01 -0.14 – 0.12
hostRacepom.sexMale -0.14 -0.23 – -0.05
Random Effects
σ2 0.44
τ00 0.07
ICC 0.86
N site 8
N ID 231
N chc 35
Observations 8085
Marginal R2 / Conditional R2 0.021 / 0.879

two way interactions

plot_model(mod, type = "pred", terms = c("hostRace", "sex"))

fixed effects

plot_model(mod, bpe = "mean", bpe.style = "dot", prob.inner = 0.4, prob.outer = 0.8)

hypothesis testing?

Manova for pomonella

summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")
##                    Df   Wilks approx F num Df den Df    Pr(>F)    
## hostRace            1 0.82295   3.3194      7    108  0.003071 ** 
## site                1 0.85799   2.5537      7    108  0.017937 *  
## sex                 1 0.52503  13.9576      7    108 8.859e-13 ***
## hostRace:site       1 0.69396   6.8041      7    108 1.085e-06 ***
## hostRace:sex        1 0.95188   0.7799      7    108  0.605479    
## site:sex            1 0.93727   1.0326      7    108  0.412735    
## hostRace:site:sex   1 0.90238   1.6691      7    108  0.124190    
## Residuals         114                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Pomonella

##           Reference
## Prediction Apple Haw
##      Apple    25  12
##      Haw      14  71
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##    0.786885246    0.503288443    0.703538181    0.855813669    0.680327869 
## AccuracyPValue  McnemarPValue 
##    0.006229331    0.844519267

LDA predicting host and sex for pomonella

##               Reference
## Prediction     Apple_Female Apple_Male Haw_Female Haw_Male
##   Apple_Female           23          2          9        1
##   Apple_Male              1          1          1        2
##   Haw_Female              8          0         32        3
##   Haw_Male                1          3          5       30
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##   7.049180e-01   5.725131e-01   6.155825e-01   7.840219e-01   3.852459e-01 
## AccuracyPValue  McnemarPValue 
##   8.606186e-13   9.110307e-01

MANOVA for Zepheria

# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df Pr(>F)
## sex        1 0.74708  0.67708      5     10 0.6507
## site       1 0.83749  0.38808      5     10 0.8461
## Residuals 14

LDA for Zepheria

##                  Reference
## Prediction        Zepheria_Female Zepheria_Male
##   Zepheria_Female               5             3
##   Zepheria_Male                 4             5
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##      0.5882353      0.1793103      0.3292472      0.8155630      0.5294118 
## AccuracyPValue  McnemarPValue 
##      0.4062810      1.0000000

LDA for sex +site zepheria

##                              Reference
## Prediction                    Zepheria_Female_EastLansing
##   Zepheria_Female_EastLansing                           2
##   Zepheria_Female_MtPleasant                            0
##   Zepheria_Male_EastLansing                             4
##                              Reference
## Prediction                    Zepheria_Female_MtPleasant
##   Zepheria_Female_EastLansing                          2
##   Zepheria_Female_MtPleasant                           1
##   Zepheria_Male_EastLansing                            0
##                              Reference
## Prediction                    Zepheria_Male_EastLansing
##   Zepheria_Female_EastLansing                         1
##   Zepheria_Female_MtPleasant                          5
##   Zepheria_Male_EastLansing                           2
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.29411765    -0.05699482     0.10313551     0.55958272     0.47058824 
## AccuracyPValue  McnemarPValue 
##     0.95792652     0.03207164

MANOVA for Mendax

summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df    Pr(>F)    
## sex        1 0.57727   6.1024      6     50 7.541e-05 ***
## site       2 0.10716  17.1232     12    100 < 2.2e-16 ***
## sex:site   2 0.74705   1.3082     12    100    0.2258    
## Residuals 55                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Mendax

##                Reference
## Prediction      Mendax_Female Mendax_Male
##   Mendax_Female            20          15
##   Mendax_Male              13          13
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.54098361     0.07072905     0.40849889     0.66935590     0.54098361 
## AccuracyPValue  McnemarPValue 
##     0.55241652     0.85010674

LDA for sex +site mendax

##                         Reference
## Prediction               Mendax_Female_Fenville Mendax_Female_OtisLake
##   Mendax_Female_Fenville                      7                      1
##   Mendax_Female_OtisLake                      1                      5
##   Mendax_Female_Sewanee                       0                      2
##   Mendax_Male_Fenville                        6                      0
##   Mendax_Male_OtisLake                        0                      1
##   Mendax_Male_Sewanee                         0                      0
##                         Reference
## Prediction               Mendax_Female_Sewanee Mendax_Male_Fenville
##   Mendax_Female_Fenville                     0                    5
##   Mendax_Female_OtisLake                     1                    0
##   Mendax_Female_Sewanee                      4                    1
##   Mendax_Male_Fenville                       0                    5
##   Mendax_Male_OtisLake                       1                    1
##   Mendax_Male_Sewanee                        4                    0
##                         Reference
## Prediction               Mendax_Male_OtisLake Mendax_Male_Sewanee
##   Mendax_Female_Fenville                    0                   0
##   Mendax_Female_OtisLake                    1                   2
##   Mendax_Female_Sewanee                     2                   3
##   Mendax_Male_Fenville                      0                   1
##   Mendax_Male_OtisLake                      1                   1
##   Mendax_Male_Sewanee                       2                   3
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##    0.409836066    0.283523654    0.285504382    0.543223627    0.229508197 
## AccuracyPValue  McnemarPValue 
##    0.001280557            NaN

MANOVA for cingulata

summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df    Pr(>F)    
## sex        1 0.78806   1.2774      4     19 0.3136059    
## site       1 0.38209   7.6817      4     19 0.0007371 ***
## sex:site   1 0.70133   2.0229      4     19 0.1319094    
## Residuals 22                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LDA for Cingulata

##                   Reference
## Prediction         Cingulata_Female Cingulata_Male
##   Cingulata_Female                5              7
##   Cingulata_Male                  7              7
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##     0.46153846    -0.08333333     0.26587122     0.66629178     0.53846154 
## AccuracyPValue  McnemarPValue 
##     0.83731354     1.00000000

LDA for sex +site cingulata

##                             Reference
## Prediction                   Cingulata_Female_SouthBend Cingulata_Female_Urbana
##   Cingulata_Female_SouthBend                          1                       2
##   Cingulata_Female_Urbana                             2                       0
##   Cingulata_Male_SouthBend                            4                       0
##   Cingulata_Male_Urbana                               3                       0
##                             Reference
## Prediction                   Cingulata_Male_SouthBend Cingulata_Male_Urbana
##   Cingulata_Female_SouthBend                        5                     2
##   Cingulata_Female_Urbana                           4                     1
##   Cingulata_Male_SouthBend                          0                     0
##   Cingulata_Male_Urbana                             1                     1
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##    0.076923077   -0.243027888    0.009455391    0.251302917    0.384615385 
## AccuracyPValue  McnemarPValue 
##    0.999943167    0.389256527

MANOVA for Cornivora

There is only one site

# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")
##           Df   Wilks approx F num Df den Df Pr(>F)
## sex        1 0.54542  0.83346      2      2 0.5454
## Residuals  3

LDA for cornivora

Cant show a plot

##                   Reference
## Prediction         Cornivora_Female Cornivora_Male
##   Cornivora_Female                2              1
##   Cornivora_Male                  0              2
##       Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull 
##      0.8000000      0.6153846      0.2835821      0.9949492      0.6000000 
## AccuracyPValue  McnemarPValue 
##      0.3369600      1.0000000